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The study of the interplay between the structure and dynamics of complex multilevel systems is a pressing 
challenge nowadays. In this paper, we use a semi-annealed approximation to study the stability properties of 
Random Boolean Networks in multiplex (multi-layered) graphs. Our main finding is that the multilevel structure 
provides a mechanism for the stabilization of the dynamics of the whole system even when individual layers 
work on the chaotic regime, therefore identifying new ways of feedback between the structure and the dynamics 
of these systems. Our results point out the need for a conceptual transition from the physics of single layered 
networks to the physics of multiplex networks. Finally, the fact that the coupling modifies the phase diagram and 
the critical conditions of the isolated layers suggests that interdependency can be used as a control mechanism. 

PACS numbers: 89.75.Fb, 89.75.Hc, 89.75.-k, 05.70.Fh 



Nearly four decades ago, Random Boolean Networks 
(RBNs) were introduced as a way to describe the dynamics of 
biochemical networks 01-01- RBNs JHHIl consider that each 
gene of a genetic regulatory network is a node of a directed 
graph, the direction corresponding to the effect of one gene 
on the expression of another. The nodes can be in one of two 
states: they are either on (1) or off (0) - i.e. in the case of a 
gene its target protein is expressed or not. The system so com- 
posed evolves at discrete time steps. At each time step nodes 
are updated according to a boolean rule assigned to each node 
that is a function of its inputs. Notwithstanding the high sim- 
plicity of RBNs models, they can capture the behavior of some 
real regulatory networks |]91 allowing for the study of several 
dynamical features, above all their critical properties. How- 
ever, although some coupled Boolean networks have been in- 
vestigated | Hi 11]. the vast majority of works has considered 
RBNs as simplex networks, in which a single graph is enough 
to represent all the interactions a given gene is involved in. 

The previous description implicitly assumes that all bio- 
chemical signals are equivalent and then collapses informa- 
tion from different pathways. Actually, in cellular biochem- 
ical networks, many different signaling channels do work in 
parallel II 1211 . i.e., the same gene or biochemical specie can be 
involved in a regulatory interaction, in a metabolic reaction 
or in another signaling pathway. Here, we introduce a more 
accurate set up for the topology of biochemical networks by 
considering that different operational levels (pathways) are in- 
terconnected layers of interaction. In terms of graphs, this 
topology is more consistent with a multiplex network lfl3l[l4ll 
(see Fig. [TJ in which each level would represent the differ- 
ent signaling pathways or channels the element participates 
in. On the other hand, accounting for the multilevel nature 
of the system dynamics also represents a point of interest by 
itself, as this allows to inspect what are the consequences of 
new ways of interdependency between the structure and the 
dynamics. In this sense, the dynamics we inspect is general 
enough so as to serve as a null model for many other complex 
dynamical processes. 



In this paper, we study the stability of Boolean networks de- 
fined at multiple topological layers. In particular, we inspect 
a Boolean multiplex network model, in which each node par- 
ticipates in one or more layers of interactions, being its state 
in a layer constrained by its own state in another layer. There- 
fore, we focus on the case of canalizing rules. Boolean func- 
tions are canalizing if whenever the canalizing variable takes 
a given value, the canalizing one, the function always yields 
the same output. Capitalizing on a semi-annealed approxi- 
mation, we analytically and numerically study the conditions 
defining the stability of the aforementioned system. By doing 
so, we show that the interdependency between the layers can 
be enough to either stabilize the different levels or the whole 
system. Remarkably, this also happens for parameter values 
where the sub-systems, if isolated, were unstable. 

Let us first define in mathematical terms the structure of the 
multiplex network of N nodes per layer and M layers in Fig. 
Q] which can be fully encoded in two objects lfl5ll . First, we 
have the NxM incidence matrix Bi a , whose elements are 1 if 
node i appears in layer a and otherwise. Secondly, we intro- 
duce an adjacency tensor, Aij a , whose elements are 1 if there 
is a link between nodes i and j in layer a and otherwise. 
With these two basic objects, one can generalize the different 
descriptors used in simplex networks. For instance, the total 
degree of node i will be K'i = £\ Q Aja = J2 a where 
Ki a is the degree of node i in layer a. Moreover, the mul- 
tilevel structure gives rise to new topological metrics that are 
not defined in single-layered networks. We define the multi- 
plexity degree of a node as the number of layers in which it 
appears as m = ^2 Bi a . Note that the number of different 
nodes in the multiplex will be then N = NAI — J2i( K i ~ -0- 

Next, let us consider a state vector 

x(t) = (xi(t),...,Xff {t) ), (1) 

where Xi(t) G {0, 1} and a set of update functions such that 

£i(t) = /i(s J -ei*( i )(t-l)). (2) 
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its neighbors in layer a and also on the state of its neighbors 
in f3 via the auxiliary function f@. Suppose that the canal- 
izing state in a and /3 is 1 (the discussion for would be 
identical). Then, the updating function of i can be written 
as fi{fi,fi) = f a V f fj , being V the Boolean operator OR. 
From the definition of the activities and the previous relation, 

it follows that E[aj] = 2~ (Kl-1) 2p(l - p), which is differ- 
ent from the value one would obtain in the case of a simple 
canalizing function. Similarly, for the sensitivity one gets 
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FIG. 1 : (color online) The multiplex network is built up by randomly 
connecting N nodes per layer. With probability a, each of the N 
nodes can be present in both layers. Therefore, the total number of 
different nodes in the system is N = (2 — a)N. In the example of 
the figure, the whole system is made up of N = 13 nodes, of which 
3 are present in the two layers and there are 5 additional nodes per 
layer, therefore N = 8 and a = 3/8. 



where T'^(i) refers to all the incoming neighbors j of node i 
at each layer a, with a = 1 . . . M, 

Equations (Q]|2| define a Boolean multilevel (or multiplex) 
graph. In addition, due to the multiplex nature of the network, 
we also define a set of update functions for each layer as 



4(t) 
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where now the arguments of the function are restricted to the 
specific layer a = I. Equation (O governs how each node is 
updated in each layer. So, Eq. (O can be rewritten as 



(4) 



where /j is a canalizing function of its inputs. These defi- 
nitions allow investigating how the stability of the Boolean 
model is affected by the multilevel structure of the system and 
by the existence of nodes with different multiplexity degrees. 

We first inspect the dependency of the average sensitiv- 
ity , which has been shown to be a useful order param- 
eter in RBNs |3 12], on the multiplexity degree m. Fol- 
lowing |16j, we write the activity a? of the variable Xj in a 
function / of K inputs as o{ = i J2> 



~j 2« ^xG{0,l} ~d^~' Where 

= /( x (j,o)) © /( x (j,i)) and X (j,-R) represents a ran- 
dom vector x G 0,1 with the jth input fixed to R and 
© is the arithmetic addition modulo 2. Similarly, assum- 
ing that the inputs are also uniformly distributed, the av- 
erage sensitivity is equal to the sum of the activities, i.e., 



* f = E*i£[x[/(x© ei )^/(x)]] = £* : 



is a zeroes vector with 1 in the i-th position, and x\M i s an 
indicator function that is equal to 1 if and only if A is true. 

To illustrate how multiplexity affects the sensitivity of a 
node, without loss of generality, we study analytically and nu- 
merically a multiplex network of two layers. Let us denote 
by p the bias of the Boolean functions, and a and j3 the two 
respective layers. Due to the multilevel nature of the interac- 
tion network, a node i in our model depends on the state of 



where E[s^ } — 2p(l — p)K a is the expected average sensi- 
tivity of a function in layer a if it were isolated. 

Next, we study the stability of the Boolean multiplex sys- 
tem using a semi-annealed approximation IU8I1 . This approach 
considers the network as a static topological object while the 
update functions f\ (I = a, (3) are assigned randomly at each 
time step. Thus, we can write the update function for the com- 
ponents of the difference vector y(f) = (| x(t) — x(t) |), 
where x is a perturbed replica of x in which a (small) fraction 
of the nodes were flipped, yielding 



W(*)=?i[l- 1(1 -%(*-!))] 



(6) 
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which is equivalent to the expression derived in 111811 . but also 
taking into account Eq. (O, with qi = 2p{l — p) for a sim- 
plex graph and Fj being the set of all neighbors of i in all 
layers. Considering a small perturbation, linearization of Eq. 
© around the fixed point solution y(t) =0 leads to 
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Vi(t + 1) 
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that can be written in matrix form as y(t + 1) = Qay(t), 
with Qij a = 2~ < - Ki ~ 1 )qiAij a . The largest eigenvalue, Aq, of 
the matrix Q = J2 a Q<* g° verns the stability of the system 
ifljll . It is worth noticing that the latter refers to the stability 
condition for the whole system and, given a fixed topology for 
each layer, it depends on the multiplexity degree Ill9ll . For 
the case of nonuniform kj we obtain an analogous mean-field 
approximation to Aq in 111 811 . 
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(2-( K >- v )q i KfKf xt ) 
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where (K) is the average degree of the multiplex. Note that 
the stability of the multiplex depends on m and K™K° ul , 
which, in general, are not independent variables — thus, q~i 
and Ki are anticorrelated. To find the critical condition let 
P(k = n) be the probability that a node in the whole sys- 
tem has multiplexity degree n. This magnitude depends on 
the same quantity but at the single layer level as P(n = n) = 
Jj^P{k — n), where P(k = n) is the probability that a ran- 
domly chosen node of a layer has multiplexity degree n. For 
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FIG. 2: (color online) Color-coded average Hamming distance for 
the whole system with fixed observed connectivity (K ) = 2.9 for 
different values of the hidden connectivity {Kh}, and the probabil- 
ity for a node in a layer to be present also in the other layer a. The 
network is composed of N = 10 3 nodes per layer as explained in 
FigE The continuos line is the solution (zeros) of Eq. Jilt . Simu- 
lations were performed for an initial Hamming distance of 0.01 and 
the results are averages over 50 realizations of the network and 300 
random initial conditions. 



the average degree of the multiplex we have: 



<*> = £ 



(M\ 



P(« = n)£<tf,) = 5r £<#,), (9) 
i i 



where (Ki) is the average degree of layer I. 

Inserting the previous expression into Eq. (|8) and consider- 
ing the case in which there are no correlations between K m 
and K out , one gets, 

i 



M -1 



with h = l...,M,h = l...,M and (q) = J2r!=i ?(« = 
n)P(n = n) is the average sensitivity on a layer. It is worth 
noticing that the first term on the l.h.s. of Eq. ([Tol l is the ex- 
pression one would obtain using an annealed approximation. 
The second term is always positive. Therefore, it captures 
the stabilizing effects of multiplexity, rightly predicting or- 
dered behavior in regions in which the annealed approxima- 
tion would not. 

Once we have derived the critical condition for a system 
made up of an arbitrary number of layers, let us compare the 
analytical results with numerical simulations for a two-layers 
system with qi = q. Let a be the probability for a node in 
a layer to be present also in the other layer, then we have 
P(k = 2) = er and P(k = 1) = 1 — a. Besides, for the 
sake of simplicity, consider that the average connectivity of 
one layer is observed, (K ), and fixed (for instance, because 
one measures it), and that the average connectivity of the other 



FIG. 3: (color online) The lines are the solution (zeros) of Eq. Jilt 
for different values of the hidden connectivity (A"h)> the observed 
connectivity (K ) and the probability of a node belongs to both lay- 
ers a. We have set qi — q — |. 



layer is unknown or hidden {Kb). Recalling that the size of 
the multiplex system is N = (2 — a)N —where N is the 
number of nodes per layer—, the mean connectivity (K) can 



be written as (K) 



(K h ) + (K ) 



, which leads to the following 



expression for the critical condition of the two-layers system 



.((K h ) + (K ))-(l-a) 



(K b ) + (K a ) 



= i (ID 

2q 



that as a function of a and (K^) gives an hyperbolic critical 
curve. 

To verify that our analytical calculations are valid, we have 
performed extensive numerical simulations of the Boolean dy- 
namics on a random multiplex network made up of two layers 
in which N nodes are randomly connected among them and 
only a fraction a of them are present on both layers. As it is 
customarily done, we test the stability of the system by mea- 
suring the long-time Hamming distance for different trajecto- 
ries generated from two close initial states. Figure|2]shows the 
results obtained when the mean connectivity (K ) of a layer is 
fixed and both a and the mean connectivity of the other layer 
(Kb) change (the Hamming distance is color coded as indi- 
cated). First, we note that the transition from stability to an 
unstable regime nicely agrees with the theoretical prediction. 
Secondly, it is worth highlighting a new effect linked to the 
multi-level nature of the system: the region of low (K^) and 
low a is unstable despite the fact that those values of (Kb) 
would make the hidden layer, in a simplex graph description, 
stable. However, due to the low coupling (er), the instability 
of the multiplex is determined by that of the observed layer, 
the leading one. Admittedly, when increasing the coupling a 
the stable (hidden) layer is able to stabilize the whole system. 

We have further explored the dependency between the sta- 
bility of the multiplex and the average degrees of both lay- 
ers. Figure [3] shows the analytical solution of Eq. (Q~T) for 
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FIG. 4: (color online) Critical curves for a network made up of 10 4 
nodes per layer as a function of the probability of a node to be part of 
both layers a, and the hidden connectivity (A'j,). The blue line corre- 
sponds to the critical curve when a single layer is observed while the 
red one refers to the whole system. The rest of simulation parameters 
are the same as for the other figures. 



Summing up, we have studied the effect of multiplexity on 
the stability of Boolean multilevel networks. In particular, we 
have addressed two important (and complementary) cases: the 
stability of the system as a whole and that of an observed 
layer which is coupled to other hidden layers. Our main re- 
sult shows that there is a region of parameters for which either 
a single layer or the whole system can be stabilized by the 
presence of another stable sub-system (layer). On more gen- 
eral grounds, the latter mechanism supports the need to study 
complex interdependent systems explicitly incorporating their 
multilevel nature. As we have shown, unexpected results can 
emerge as a consequence of new ways of feedback between 
the structure and the dynamics of such systems, including the 
possibility of using interdependency to control the stability of 
a system. 
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different values of (K^) and (K ). The results show a very 
rich phase diagram. Depending on the values of both con- 
nectivities, a double transition from a chaotic regime to an 
ordered one and again to another chaotic regime is predicted. 
More interestingly, the transition from the ordered to the dis- 
ordered regime does not depend on a only when both lay- 
ers operate at their respective critical points, namely, when 
(K h ) = (K ) = l/q = 2. 

Up to now, we have analyzed the stability of the multiplex 
system. In practice, it is more common to have access to 
only one layer, so that one can measure the stability of that 
layer given that it is connected to a hidden (inaccessible) one. 
Therefore, it is also important to inspect the stability condi- 
tion of a single layer within the multiplex. To this end, we 
should solve Eq. dS) taking into account only the nodes that 
belong to the layer whose stability is scrutinized. In this case, 
the critical condition reads 

\ ({K h f - (K f + 2(A h ) { K o) ) + = <*■>+»<*»> . 

(12) 

Figure [4] compares results of simulations for a larger network 
of N = 10 4 nodes per layer with the theoretical solution (Eq. 
( fT2l . blue line) showing again a good agreement between an- 
alytical and simulation results. Remarkably, the results show 
that a single ingredient —the multilevel nature of the system 
— can explain why there are biologically stable systems that 
are however theoretically expected to operate in the unstable 
regime (i.e., their average degree is larger than 1 / q). In other 
words, the sole reason could be that these systems are not iso- 
lated, but are coupled to other hidden layers that, if ordered, 
can stabilize the system. Finally, for the sake of comparisons, 
we have also represented in Fig.|4](red line) the case shown in 
Fig.|2]but for the same larger system size. 
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